#!/usr/bin/perl -w

$fin1 = "Cavity_Cluster.txt";
$fin2 = "all.one.line.fa";
$fin3 = "aa_prop.txt";
$fin4 = "epi_sites.txt";
$fout = "predSet.txt";

$main_strain = "ABQ44394"; # New Jersey...

open FIN1, $fin1;
open FIN2, $fin2;
open FIN3, $fin3;
open FIN4, $fin4;
open FOUT, ">", $fout;

$epi_name = ["Sa", "Sb", "Ca", "Cb", "Pa", "Pb"];

# epitope: Sa, Sb, Ca, Cb, Pa, Pb
$dataset = {
	"Ndon" => [0, 0, 0, 0, 0, 0],
	"Nacc" => [0, 0, 0, 0, 0, 0],
	"Npc" => [0, 0, 0, 0, 0, 0],
	"Nnc" => [0, 0, 0, 0, 0, 0],
	"Nh" => [],
};

=start
$dataset = {
	"Ndon" => [0, 0, 0, 0, 0, 0],
	"Nacc" => [0, 0, 0, 0, 0, 0],
	"Npc" => [0, 0, 0, 0, 0, 0],
	"Nnc" => [0, 0, 0, 0, 0, 0],

	"Nh1" => [0, 0, 0, 0, 0, 0],
	"Nh2" => [0, 0, 0, 0, 0, 0],
	"Nh3" => [0, 0, 0, 0, 0, 0],
	"Nh4" => [0, 0, 0, 0, 0, 0],
	"Nh5" => [0, 0, 0, 0, 0, 0],
	"Nh6" => [0, 0, 0, 0, 0, 0],
};
=cut

# TEST ref length
#$tt = [1, 2, 3];
#print @$tt + 0; exit; # get length

# 20 amino acids hash
$aa_prop = {
	# "A" => [],
	# "R" => 
};

$seqs = {}; # store aa seqs
# $train_set = {}; # store pair diff & distance  ==> write to file only.
$epi_sites = {};
$cavity = {};

$tmp = <FIN3>; # drop the title
while(<FIN3>)
{
	# print $_;
	# exit;#
	chomp;
	
	/(\w)\t(\d)\t(\d)\t(\d)\t(\d)\t(.*)/;
	#print $1, $2, $3, $4, $5, $6;
	# Res, don, acc, pc, nc, h
	$aa_prop->{$1} = [$2 + 0, $3 + 0, $4 + 0, $5 + 0, $6 + 0];
}

# TO DO with "-"  ===> GAP
$aa_prop->{"-"} = [0, 0, 0, 0, 0]; ###??

print $aa_prop->{"A"}->[4];
print $aa_prop->{"N"}->[4];
#exit; ##

# store seqs
while(<FIN2>)
{
	chomp;
	# substr... # delete '>'
	/>(.*)/;  # BETTER than 'substr'
	#print $1; exit;
	$tmp = <FIN2>; 
	chomp $tmp;
	$seqs->{$1} = $tmp; # strain name => seqs
	
	#print $seqs->{$1}; exit;
}

# epi_sites
$tmp = <FIN4>;
while(<FIN4>)
{
	chomp;
	/(\w*)\t(\d*)\t(.*)/;
	$epi_sites->{$1} = [$2 + 0, $3]; # num of sites, every site
}
#print $epi_sites->{"Sa"}->[0]; exit;

# cavity
$cav_count = 0;
$epi_n = "";
my @sites;
while(<FIN1>)
{
	chomp;
	$tmp = $_;
	if($tmp =~ />/)
	{
		$cav_count = 0; # IMPORTANT
		#print; exit;
		$tmp =~ />(.*)/;
		$epi_n = $1;
		print $1; #exit;
		#next;
	}
	else
	{
		#print $tmp, "\n";
		$tmp =~ /cavity: (.*)/;
		print $1, "\n"; #exit;
		# @sites = split /\t/, $1; # it's space, not '\t'
		@sites = split / /, $1;
		# print $sites[2], "\n";
		
		
		#print $epi_name, "\n";
		# $cavity->{$epi_name}->[$cav_count] = \@sites; # REFERENCE 
		# ==> PROBLEM: all the address is the same
		
		# ANONYM ==> NO problem
		$cavity->{$epi_n}->[$cav_count] = [@sites]; # ANONYM
		#print $cav_count, "\n";
		print $sites[1], "\n";
		print $cavity->{$epi_n}->[0]->[1], "\n";
		
		$cav_count++; # IMPORTANT !
	}
}
#print defined($cavity->{"Sa"}); exit; ##
#print defined($cavity->{"Sa"}->[6]); exit; ##
#print defined($cavity->{"Sa"}->[5]); exit; ## ==> use to know how much cavities .... in a epitope
print "length: ", @{$cavity->{"Sa"}} + 0, "\n"; # 6

#print FOUT "Ndoc\tNacc\tNpc\tNnc\tNh\t\n";

# construct the data set
foreach $strain2 (keys(%$seqs))
{
	$strain1 = $main_strain;
	print "\n", $strain1, "\t", $strain2; #exit; ##
	
	$ds = get_pair_data($strain1, $strain2);
	
	# print $ds->{Ndon}->[0]; # the SAME
	#print $ds->{"Ndon"}->[0], "\t";
	print "\t", $ds->{"Nh"}->[0]->[0], "\t";
	
	# write => Method 1
	for($i =0; $i < 6; $i++)
	{
		# ALL
		print FOUT $ds->{"Ndon"}->[$i], "\t", $ds->{"Nacc"}->[$i], "\t", 
			$ds->{"Npc"}->[$i], "\t", $ds->{"Nnc"}->[$i], "\t";
		$len = @{$ds->{"Nh"}->[$i]} + 0;
		$tmp = $ds->{"Nh"}->[$i];
		
		for($j = 0; $j < $len; $j++)
		{
			print FOUT $tmp->[$j], "\t";
		}
		
		
		# Sa
=start
		print FOUT $ds->{"Ndon"}->[0], "\t", $ds->{"Nacc"}->[0], "\t", 
		$ds->{"Npc"}->[0], "\t", $ds->{"Nnc"}->[0], "\t";
		$len = @{$ds->{"Nh"}->[0]} + 0;
		$tmp = $ds->{"Nh"}->[0];
		
		for($j = 0; $j < $len; $j++)
		{
			print FOUT $tmp->[$j], "\t";
		}
=cut	
		# Sb
		
		# Ca
		
		# Cb
		
		# Pa
		
		# Pb
		
	}
	print FOUT $strain1, "\t", $strain2; ##
	
	# Method 2=>
		
	# new line
	print FOUT "\n";
}

sub get_pair_data
{
	my($strain1, $strain2) = @_;
	# my($Ndon, $Nacc, $Npc, $Nnc, $Nh); # INIT them every loop
	# find the seqs, respectively.
	
	# if seqs in the file 'combineSeq.fa' ...
	#print $strain1, "\n";
	#print $seqs->{$strain1};
	#exit;
	
	@str1 = split(//, $seqs->{$strain1});
	@str2 = split(//, $seqs->{$strain2});
	
	#print "\n", $strain1, "\n", $strain2;
	#print "\nstr1=", @str1, "\nstr2=", @str2; exit;
	
	for($i = 0; $i < 6; $i++)
	{
		
		my($Ndon, $Nacc, $Npc, $Nnc, $Nh) = (0, 0, 0, 0, 0);
		
		#print $i; exit;
		#print $epi_name; exit;
		#print $epi_name->[$i]; exit;
		$tmp = $epi_sites->{$epi_name->[$i]};
		#$num = $tmp->[0];
		@sites = split(/ /, $tmp->[1]);
		#print "\n", $tmp->[1], "\n", "@sites"; exit;
		
		$cav_site = $cavity->{$epi_name->[$i]};
		
		#print $str1[0]; exit;
		foreach(@sites)
		{
			#$prop1 = $aa_prop->{$str1[$_ + 0]}; <====> start from '0', not '1', must '-1'
			#$prop2 = $aa_prop->{$str2[$_ + 0]};
			#print $_, "\t", $str1[$_ + 0], "\t", $str2[$_ + 0], "\n";
			$prop1 = $aa_prop->{$str1[$_ - 1]};
			$prop2 = $aa_prop->{$str2[$_ - 1]};
			#print "\n", $_, "\t", $str1[$_ - 1], "\t", $str2[$_ - 1];
			
			#need to deal with "-"
			
			#print $prop1->[0]; #exit; #
			#print $prop2->[0]; exit; #
			
			$Ndon += $prop1->[0] - $prop2->[0];
			$Nacc += $prop1->[1] - $prop2->[1];
			$Npc += $prop1->[2] - $prop2->[2];
			$Nnc += $prop1->[3] - $prop2->[3];
		}
		#exit; ##
		$dataset->{"Ndon"}->[$i] = $Ndon;
		$dataset->{"Nacc"}->[$i] = $Nacc;
		$dataset->{"Npc"}->[$i] = $Npc;
		$dataset->{"Nnc"}->[$i] = $Nnc;
		
		$j = 0;
		print "\nlength: ", @{$cav_site->[$j]} + 0;
		while(defined($cav_site->[$j]))
		{
			$cav = $cav_site->[$j];
			$len = @$cav + 0;
			print "\n@$cav";
			
			for($k = 0; $k < $len; $k++)
			{
				# WRONG for the same reason above: '-1'
				#$prop1 = $aa_prop->{$str1[$cav->[$k] + 0]};
				#$prop2 = $aa_prop->{$str2[$cav->[$k] + 0]};
				
				$prop1 = $aa_prop->{$str1[$cav->[$k] - 1]};
				$prop2 = $aa_prop->{$str2[$cav->[$k] - 1]};
				#print "\n", $k, "\t", $str1[$cav->[$k] - 1], "\t", $str2[$cav->[$k] - 1]; #exit;
				
				$Nh += $prop1->[4] - $prop2->[4];
			}
			
			$dataset->{"Nh"}->[$i]->[$j] = $Nh;
			print "\n", $Nh; #exit;
			
			$j++;
		}
		#exit;
	}
	#exit;
	
	return $dataset;
}


